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We investigate the nonlinear dynamics of a mesoscopic driven Duffing oscillator in a quantum regime. In 
terms of Wigner function, we identify the nature of state near the bifurcation point, and analyze the transient 
process which reveals two distinct stages of quenching and escape. The rate process in the escape stage allows 
us to extract the transition rate, which displays perfect scaling behavior with the driving distance to the bifurca- 
tion point. We numerically determine the scaling exponent, compare it with existing result, and propose open 
questions to be resolved. 
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A broad class of physical systems such as Josephson junc- 
tion, trapped electron or ion, and nano-mechanical oscillator, 
can be well described by the Duffing oscillator under proper 
conditions. One of the most profound features of a driven 
Duffing oscillator (DDO) is the dynamical bifurcation. Near 
the bifurcation point, the oscillator state is highly sensitive 
to perturbation. This property can be exploited for appli- 
cations such as sensing device, amplifier, and logic device. 
Most recently, for instance, the superconductor circuit based 
on Josephson junction has been exploited for quantum mea- 
surement of superconducting qubits Hl-fl]. This device is 
termed as Josephson bifurcation amplifier (JBA), holding ad- 
vantages such as fast speed, high sensitivity, low backaction, 
and absence of on-chip dissipation. 

Despite that the classical bifurcation of DDO is well- 
known, the quantum dynamics in the bistable region and near 
the bifurcation point has been a new and significant subject 
in the past years fs Uloll . This new trend is motivated mostly 
by the advent of approaching the quantum regime of nano- 
mechanical oscillators, as well as the bifurcation-based quan- 
tum measurement devices. For instance, the quantum signa- 
ture in the bistable region of a DDO was proposed, based 
on simulating a Lindblad-type master equation and compar- 
ing the Wigner function with classical probability distribution 
in phase space In terms of amplitude and phase responses 
to the driving frequency, quantum behaviors of DDO such as 
resonant tunneling and photon-assisted tunneling were also 
discussed |g]. Moreover, in Ref. switching rate be- 

tween the bistable states near the bifurcation point, due to 
quantum and/or thermal fluctuations, is estimated by means of 
the WKB theory or semiclassical methods such as the mean- 
first-passage-time approach. 

In this letter we consider a mesoscopic DDO, with about 
more than ten levels that is in between the quantum few-level 
and the classical dense-level (or continuum) limits. In this 
regime, the quantum effect is apparently significant. How- 
ever, at the same time, how the DDO's nonlinearity manifests 
itself is unclear and of interest, since the few-level (e.g. 2- or 



3-level) system should have no such behaviors as bistability 
and bifurcation. Our present study will demonstrate the ex- 
istence of bistable region, characterize the quantum nature of 
the states, and investigate the quantum transition near the clas- 
sical bifurcation point which displays perfect and new scaling 
behavior with the driving strength. 

Model and Method. — The Duffing oscillator in the pres- 
ence of driving is described by the Hamiltonian 
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H (t) = + -mtfx 2 
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yt + F(t)x. 



(1) 



For the JBA setup, F(t) = Fo(e' yt + e~ iyt ) describes the mi- 
crowave driving. Other parameters are related to the JBA cir- 
cuit quantities as: m = (h/2e) 2 C, Q. = -\/2eI c /(hC), Fq = 
hl/(2e),and y = ra£2 2 /24; with C the capacitance of the 
Josephson junction, I c the critical current, and / the driving 
current. In this context, x denotes the phase difference across 
the Josephson junction. 

In addition, the Duffing oscillator is affected by environ- 
ment, which together with the coupling can be modelled as 
He = Zi[ J w i a»?^i 2 /2 + /) 2 /2m,] - + * 2 Z; A 2 l(2m i u 2 ). 
Typically, the spectral density of the bath, J(a>) = n A 2 S{lo- 
w;)/(2m,Wi), in Ohmic case reads J(u>) = niKcaexp(-a> / a> c ), 
with k the friction coefficient, and ai c the high frequency cut- 
off. For later use, we also introduce b = Yii^&il V2, with 
hi — (tnjLOiXj + ipi) I y2rn^ha>i. 

In the weak coupling limit to the environment and un- 
der Markovian approximation, the dissipative dynamics of 
the DDO is governed by the quantum master equation (see 
Ref. Ill ill for more details) 



p(t) = - l -[H(t\p(t)] - i{[f, Qp(t)] + H.c. 
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Here, p(t) is the reduced density matrix of the oscillator; 
H(t) = H (t) + x 2 mKO) c /7T, and Q = [C(-£) + C(-£)]x/2. 
The Liouvillian £. is defined through its action on an arbitrary 
operator O as: £.0 = hT l [H{t) - F{t)x,0]. The superop- 
erators C(£.) and C{£.) are the Fourier transform of the bath 
correlation functions: C(_£) = J ^ dtC(t)e lJ - t , and C{£) = 
J m dtC(i)e tU . The correlators C(t) and C(t) are defined by 
C(t) = Tr £ [^ t (f)^(0)p £ ], and C(t) = Tr E [b(t)bH0)p E ], where 
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Pe is the thermal-equilibrium density operator of the environ- 
ment. 

Qualitative Considerations. — In the absence of driving, the 
Duffing oscillator described by Eq. (Q~|i has only finite number 
of bound states. This can be seen from the potential profile, 
V(x) = m£l 2 x 2 /2 - yx 4 , which defines a single well with iden- 
tical barrier height V = ra 2 Cl 4 /(16y) at x = ± \Jm£l 2 /(4y). 
As a rough estimate, the number of bound states is the ratio of 
V and Ml, which gives N = m 2 Q 4 /(16ySQ) = mQ/(16Sy) = 
N/(16y). We will see later that S = m£l/h defined here is 
an important characteristic quantity. We also introduced a 
reduced nonlinear coefficient, y = y/(mQ. 2 ). In our model, 
y = raQ 2 /24, so approximately the number of bound state is 
3S/2. In the experiment of Ref. Jit], S - 366, which implies a 
classical DDO. In the present work, we consider a mesoscopic 
regime, by assuming possible parameters I c = 39nA, C 
0.9 1 pK K = O.Olfl, and cj c = 10Q. Accordingly, S = 12. 

Under proper conditions J^t], the DDO exhibits the most 
profound phenomenon known as bifurcation. To determine 
the bifurcation point, we present an analysis in the rotating 
phase space. Starting with the Hamiltonian Ho(t), we intro- 
duce a unitary transformation, U = exp[/vf(a"'"a)L where a is 
the annihilation operator of the Duffing oscillator. Under the 
rotating wave approximation (RWA), we obtain Jgt] 
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where 6=1- v/Q, in = m/5, and Q = Q.6. In phase space, 
the extremal points of satisfies 



0, 
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For this cubic equation, the discriminant reads A = 
[-2F /(3y x 2)] 2 + [-2mQ 2 /(3y x 3)] 3 . If and only if A < 0, 
there exist three different real roots. This implies 



2 2m 3 Q 6 ,„ 8 V3mf2 2 (5 3/2 
F Q < -( ) 1/2 = = F c . 
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This result, based on the existence condition of multiple 
steady states, coincides with that from the singularity anal- 
ysis [3J. The advantage of the present method is its appli- 
cability to more general situation, e.g., the mesoscopic case 
under present study, in which we will see that the singular bi- 
furcation is absent. Nevertheless, in this work we still refer to 
the estimated F c of Eq. ® as bifurcation point for description 
convenience, particularly as a reference value for the driving 
strength. 

Moreover, still classically, to realize the bifurcation it was 
found in Ref. [3] that the driving frequency v should be lower 
than Q and satisfy 5 > V3/(20, where Q = 0//c. Based on 
the above analysis, here we can further set a upper limit to the 
detuning. For Fo < F c , Eq. has three different real roots, 
with the largest one as X max = 2[F e /(3y)] 1 ^ 3 cos6'/3, where 
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FIG. 1: (a) and (b): Bistable feature visualized by the transient be- 
havior of x{t) =Tr[xp(t)]. Parameters: driving strength Fo = 0.SF c , 
and driving frequency v = 0.94Q. (c): Two successive stages (i.e. 
quenching and escape) towards the LAS, for driving near the critical 
point as exemplified here by Fo = 0.95F r . 
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FIG. 2: (a) Phase diagram of the oscillation amplitude of stable state 
against the driving strength. The red and black curves are for initial 
states of Gaussian wavepackets centered at x = and x = 1, respec- 
tively, (b) Classical counterpart of (a), showing sharp bifurcation 
behavior. 



8 = arctan -<J(F c /Fq) 2 - 1. Obviously, this largest amplitude 
should not overcome the potential barrier of the Duffing oscil- 
lator. This consideration leads to the following inequality: 
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In our model, y = mQr/24. So we get 5 < 18/64 « 0.28. 
This is not the optimized value. More accurate consideration 
can result in even smaller upper limit. 

Quantum Dynamics near the Bifurcation Point. — Based on 
a direct simulation of the master equation, we show in Fig. 
1 the evolution of x(t) = Tr[xp(t)]. First, in Fig. 1(a) and 
(b), we demonstrate the bistable nature, for driving strength 
Fo = 0.8F C as an example. We consider two initial condi- 
tions: the ground state, and a coherent state centered at x — 1 . 
Indeed, for the mesoscopic DDO, here we find quantum me- 
chanically that the steady state does exhibit bistable behav- 
ior, say, depending on the initial condition, it arrives at ei- 
ther a small amplitude state (SAS), or a large amplitude state 
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(LAS). Nevertheless, as we will understand later in the follow- 
ing study, the results in Fig. 1(a) and (b) are not the fundamen- 
tal SAS and LAS, but their mixture with different population 
probabilities depending on the initial conditions. For driving 
not very close to F c , the "steady state" population is formed 
in relatively short time, and later the fluctuation-induced tran- 
sitions between the fundamental SAS and LAS are negligibly 
weak. 

In contrast, as shown in Fig. 1(c), for driving closer to F c 
(e.g. Fq = 0.95F C ) we find that the entire process contains 
two distinct stages, say, a (fast) quenching stage, and a suc- 
cessive (slow) escape stage. In the quenching stage, the os- 
cillator rapidly arrives at the SAS. Then, it is followed by a 
rate process (transition) to the LAS, which may be termed as 
the Kramers escape process 112I1 . Conventionally the escape 
is caused by thermal fluctuations, as described by for instance 
the mean-first-passage-time approach or Fokker-Planck equa- 
tion InRef. J2O, quantum-fluctuation induced transition 
was also investigated, by using the WKB method. The advan- 
tage of our present numerical simulation allows to formulate 
a way to extract the transition rate under more general condi- 
tions, say, in the presence of both thermal and quantum fluc- 
tuations, going beyond the existing results in limiting cases. 
This will be detailed in latter analysis. 

In Fig. 2(a) we extract data from numerical simulation as 
shown in Fig. 1 to plot the phase diagram, say, the oscil- 
lation amplitude of steady state against the driving strength, 
which shows the desirable hysteresis behavior. However, 
compared to its classical counterpart as schematically shown 
in Fig. 2(b), two differences should be remarked, (i) The 
singularity of transition from the bistable region to the sin- 
gle LAS or SAS disappears, although the transition point, i.e., 
F c — (8 y/3mQr6 3 / 2 )/9, is approximately preserved. The ba- 
sic reason for this gradual transition behavior is that for the 
present mesoscopic DDO, the stable state is a statistical mix- 
ture of the SAS and LAS in the absence of noise (i.e. thermal 
and quantum fluctuations), (ii) A dip appears in the red curve. 
Our numerical simulation shows that the time-dependent os- 
cillations of the SAS and LAS are out of phase (i.e. with a 
phase difference about n), and the SAS and LAS themselves 
depend on the driving strength nonlinearly. As a consequence 
of interplay of these two factors, the dip is formed as we ob- 
served. 

Let us proceed our further analysis with the help of 
the Wigner function, which is defined as: W(x,p,t) = 
l/(nh) f ^ {x + x'\p{t)\x - x') exp(-i2px' /h)dx' . The Wigner 
function is widely used in broad context of physics, with an 
intuitive interpretation of probability in the phase space. In 
Fig. 3(a) we show the Wigner function of the oscillator at time 
t = 160 * (2n/£2), for driving strength Fq = 0.9F C , and with 
the ground state as the initial condition. Time dependently, 
the Wigner function is in fact rotating with the driving fre- 
quency in phase space, along the classical trajectory but with 
additional diffusion because of the thermal and/or quantum 
fluctuations. 

In the transient process, after certain duration time to be dis- 
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FIG. 3: (a) Wigner function at t = 160 * (In/Q.), for driving 
F = 0.9F C and starting with the ground state. The two separated 
wavepackets correspond to the SAS and LAS, and the SAS is nearly 
a perfect coherent state, (b) Occupation probability of the LAS (P2) 
against the driving strength. In the region near F c , P2 shows expo- 
nential dependence of the driving strength. 



cussed below, the oscillator can be well described by a mixed 
state, formally as 

W(x, p,t)=P 1 (t)Ws (x, p, t) + P 2 (t)W L (x, p, f). (6) 

Here, Ws (x, p, f) and Wl(x, p, t) are, respectively, the Wigner 
functions of the intrinsic SAS and LAS associated with the 
given driving strength, but not the averaged ones shown in 
Fig. 2(a). In essence, the SAS and LAS are two limit cycles, or 
attractors, each with its own attraction basin fil llal . Pi(t) and 
Piit) are the occupation probabilities of the SAS and LAS. In 
Fig. 3(b) we plot Pi(0 at f = 160 * (In/Q) versus the driving 
strength, for two temperatures. 

Closer inspection indicates that Ws (x, p, t) is quantum me- 
chanically a pure state, i.e., Trp^ AS ^ 1, where p$As is the 
density matrix of the SAS. Moreover, it is nearly a perfect co- 
herent state. This result can be understood as follows. For 
a harmonic oscillator under the interplay of driving and dis- 
sipation, such as the optical cavity field under excitation and 
photon-loss, the steady state is exactly a coherent state 1 1-411 - 
Then, for the present nonlinear Duffing oscillator, since the 
SAS is not far from the oscillator origin, it is thus approxi- 
mately governed by a harmonic oscillation. For LAS, how- 
ever, which is far from the origin, nonlinearity is prominent, 
which makes Wl(x, p, t) not at all a coherent state, but a mixed 
state with partial coherence. 
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FIG. 4: Perfect scaling behavior hidden in the transition rate with the 
driving distance to the critical point, say, 17 = F\ — Fq. The triangles 
stand for data from numerical simulation, and the straight lines are 
linear fits. The result shows that R cc rf, and a = 1. Parameters: 
K = 0.01; r=5mK in (a), and T=50mK in (b). 



Transition Rate. — Now we return to the transient dynamics 
and focus on the escape stage as indicated in Fig. 1(c). This 
stage is described by a rate process: 



dPi(t) dP 2 (t) 

-K\P\ + K2P2, 



dt 



dt 



-K2P2 + K\P\ 



(7) 



where k\ is the escape rate from SAS to LAS, and K2 vice 
versa. In what follows, we formulate a way to determine the 
escape rate, by contacting Eq. (0 with the numerical simula- 
tion. First, the solution of Eq. (O reads 

P «) =^ + h(0)-^]e X p^' + ^ (8) 

p 2® =^ + N0)-^]exp-^'. (9) 

Pi(0) and P 2 (0) are certain "initial" values of the population 
probabilities in the escape stage. Rather than Pi(0) and P 2 (0), 
which are not well defined, we use the probabilities P 2 {tj) at 
three time points, and simply assume t$ - t 2 = t 2 — t\ = At. 
Then, based on Eqs. (O and (0, we obtain 



[Piih) - (K - \) 2 Piih)] lnQg-1) 
[\-{K-\f]At 



\n(K - 1) 
k 2 = K U 



(10) 
(11) 



where K = [P 2 (t 3 ) - P 2 {h)y[P 2 (t 2 ) ~ P 2 (h)l 

Below we focus on /q and formally assume K\ = Ce~ R/A . 
Here C is an irrelevant prefactor, and the exponential form of 
~ e~ R l x is associated with an effective activation process. In 
limiting cases, such as for classical thermal activation, R is 
the activation energy and A the temperature; while for quan- 
tum tunneling through a barrier, R is the tunneling action and 
A the Plank constant. Our present situation is a generaliza- 
tion, i.e., quantum-dynamical-tunneling dominated but also 
thermal-activation involved. So, we may view R as an effec- 
tive activation energy and A an effective Planck constant or 
temperature. 

Physically, we should expect that the transition rate de- 
pends on the driving distance to the critical point F c , i.e., 



T] = - Fq, since closer to the critical point, more easily 
can the transition to the LAS take place. Strikingly, Figure 4 
displays a perfect scaling behavior for this dependence. As- 
suming R oc if', our precise numerical fitting gives a 1. We 
noticed similar scaling behavior was found by Dykman @], 
but where a scaling exponent a = 3/2 was found instead. 

Since our simulation is for a mesoscopic DDO with a bit 
more than ten levels involved in the dynamics, we postulate 
that the scaling exponent a = 3/2 is not universal. As in 
and ifioll . our present simulation does not account for the driv- 
ing field in the dissipation terms. Although this kind of treat- 
ment is well accepted in vast areas, there is indeed some coun- 
terexamples, for instance, see the most recent Ref. Il5tl . Nev- 
ertheless, by transforming the system to a rotating frame to 
account for the driving in the dissipation terms and calculat- 
ing the fidelity of the SAS, we actually discovered the same 
scaling exponent llal . 

We noticed that in Ref. |H3l, scaling behavior of the tran- 
sition rate with the driving frequency (but not the driving 
strength) was analyzed to give a 1.3 ~ 1.4, by a rough 
fitting from a few experimental data. Meanwhile, in the 
experiments by Siddiqi et al. iflitl . an effective potential 

n r 9~l 3 / 2 
with a barrier height scaled as &U~. oc 1 — (Fo/F c ) \ , 

was employed to analyze their measured data by means of 

the thermal-activation rate oc exp(— AC/^ '■ /kaTj. Based on 

the same effective potential, a rough WKB analysis should 

result in a smaller scaling exponent for quantum rate oc 

exp^- ^AUj yn a/hj, with a the effective width of the barrier. 

It seems that our above result a 1 is in qualitative agreement 
with this analysis. Therefore, stronger experimental evidences 
and more rigorous theoretical investigations will be helpful to 
clarify this interesting issue, particularly with an extension to 
the mesoscopic regime as we estimated earlier in this work. 

Summary and Discussion. — We have investigated the 
quantum nonlinear dynamics of the DDO in a mesoscopic 
regime with a few more than ten energy levels involved in the 
dynamics. We demonstrated for the first time that the quantum 
nature significantly modifies the classical sharp-switching be- 
havior near the bifurcation point. In terms of Wigner function 
the state near the bifurcation point was identified to be a mix- 
ture of low- and high-amplitude states, and the low-amplitude 
component is nearly a perfect coherent state. More interest- 
ingly, near the bifurcation point the transient dynamics reveals 
distinct stages of quenching and escape. The latter is well 
characterized by a rate process, and the numerically extracted 
rate displays perfect scaling behavior with the driving distance 
to the bifurcation point. 

The quantum predictions of this work, in particular the scal- 
ing exponent that differs from the existing result raise 
an open question and deserve further investigations. The sta- 
tistically mixed nature of the state near the bifurcation point 
may also partially explain the discrepancy between the clas- 
sical prediction and the measurement data in the JBA experi- 
ment J2] • Related to this, the mesoscopic DDO may also sup- 
port quantum weak measurement in the transient stage, where 
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qubit state can be updated using a generalized Bayesian rule. 
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